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SUMMARY 


Four explicit finite-difference techniques designed to solve the time- 
dependent, compressible Navier-Stokes equations have been compared. These 
techniques are (1) MacCormack, (2) modified Du Fort-Frankel , (3) modified 
hopscotch, and (4) Brailovskaya. The comparison was made numerically by 
solving the quasi-one-dimensional Navier-Stokes equations for the flow in a 
converging-diverging nozzle. Solutions with and without standing normal 
shock waves were computed for unit Reynolds numbers (based on total condi- 
tions) ranging from 45374 to 2269. The results indicate that all four 
techniques are comparable in accuracy; however, the modified hopscotch scheme 
is two to three times faster than the Brailovskaya and MacCormack schemes and 
three to six times faster than the modified Du Fort-Frankel scheme. 


INTRODUCTION 


Recently, considerable interest has surfaced in the numerical solution 
of the compressible Navier-Stokes equations (refs. 1-4). Explicit numerical 
techniques have been used in most of these studies, especially those involving 
shock waves. The limited use of implicit methods is due to (1) coding com- 
plexity associated with the Navier-Stokes equations, (2) limited success in 
obtaining the large time steps as predicted by linear stability analysis, and 
(3) limited success in capturing shock waves. Another factor involved is the 
apparent success of explicit methods over implicit methods for adapting to 
the new fourth generation computers (STAR 100 and ILLIAC IV). 

The purpose of the present study is to investigate the relative merits 
of several explicit finite-difference techniques for solving the compressible, 
time-dependent Navier-Stokes equations. Some of the important aspects 
evaluated are (1) computational speed, (2) numerical accuracy, (3) computer 
storage requirements, (4) Reynolds number limitations, and (5) effects of 
artificial smoothing. The four numerical techniques investigated are (1) 
modified hopscotch, (2) MacCormack, (3) modified Du Fort-Frankel, and (4) 
Brailovskaya. Each of these methods has been used to solve a quasi-one- 
dimensional converging-diverging nozzle problem. Solutions with and without 
standing normal shock waves are presented for unit Reynolds numbers ranging 
from 45374 to 2269. 
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SYMBOLS 


O 

A nozzle cross-sectional area, m 

c speed of sound, m/sec 

O 

E total internal energy per unit volume, N-m/m 

M Mach number 

2 

p pressure, N/m 

R Reynolds number per unit length, m“^ 

S smoothing term 

t time, sec 

T temperature , K 

u velocity, m/sec 

X distance along nozzle axis, m 

At time increment, sec 

Ax space increment, m 

3 

p density, kg/m 

Subscripts: 
i space index 

t total conditions 

Superscript: 
n time index 


GOVERNING EQUATIONS AND TEST PROBLEM 


The converging-diverging nozzle problem used in this study represents a 
rigorous test case for the numerical techniques. The steady-state flow field 
is initially subsonic, goes sonic at the throat, passes through a standing 
normal shock wave in the diverging portion of the nozzle, and exits the 
nozzle with subsonic flow. Cases which do not contain a standing normal 
shock wave are also computed. For these cases the flow field downstream of 
the throat is supersonic. Different exit boundary conditions are required 
for each of these cases and will be discussed in a subsequent section. 
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The time-dependent, quasi -one-dimensional flow of a compressible, 
viscous fluid is governed by a set of three partial differential equations 
expressing the conservation of mass, momentum, and energy. These equations 
in conservative form are as follows: 
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The coefficients of viscosity (y) and thermal conductivity (k) are given by 
Sutherland's viscosity law and a constant Prandtl number assumption. 


Differencing Schemes 

Several characteristics are common to each of the numerical schemes 
evaluated in this study. They are all second-order-accurate (for the steady- 
state solution) finite-difference techniques which solve the time-dependent 
form of the governing equations in search of a final steady-state solution. 
The methods are explicit, and hence, easily programmed. In particular, the 
methods evaluated here have been chosen especially with regard to programming 
simplicity for the Navier-Stokes equations. 


Modified Hopscotch 

The current version of hopscotch was first introduced in reference 5 
where it was applied to the compressible Navier-Stokes equations for a shear 
layer mixing problem. This modified hopscotch technique (applied to eq. (1)) 
is expressed in two sweeps given by 

first sweep (i+n even) 
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second sweep (i+n odd) 
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Gottlieb and Gustafsson (ref. 6) have investigated the stability of the 
current version of hopscotch and have found it to be governed by the following 
CFL condition: 


^^CFL = 


Ax 


+ c 


(5) 


In addition reference 6 found the viscous stability condition for the 
modified hopscotch technique to be 


y At < , 
P (Ax)^ 


( 6 ) 


Gourlay (ref. 7) suggested a simplification to the standard two-sweep 
hopscotch scheme which almost entirely removed the first sweep where equation 
(3) is replaced by 


= 2U'J - U'?"'' (i+n even) (7) 

Numerical tests were performed with and without the use of equation (7) 
yielding identical results. The use of equation (7) increases the speed of 
the modified hopscotch technique by a factor of two without requiring 
additional storage. 


Modified Du Fort-Frankel 


The current version of the Du Fort-Frankel scheme was introduced by 
Gotti eib and Gustafsson (ref. 8) as follows: 
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where the last term is a stabilizing term. The value of to (stabilizing 
coefficient) must be determined by numerical experiment. The stabilizer takes 
the place of the time averaging appearing in the viscous terms of the standard 
Du Fort-Frankel scheme. This simplifies the resulting numerical code, 
especially for the Navier-Stokes equations in multiple dimensions. 


In addition to equation (8) an additional dissipative term must be 
added for stable operation 
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where the constant e is determined numerically. 


MacCormack 


The version of the two-step Lax-Wendroff scheme used in this study was 
first introduced by MacCormack (ref. 9). Using the MacCormack technique to 
difference equation (1) yields 
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where the overbar on the n superscript indicates a predicted value. The 
stability requirement for this scheme is the CFL condition (eq. (5)). In 
addition, a stability condition due to viscous effects is also present: 


y At < X 

p(Ax)2 “ 2 (12) 


To surpress pointwise oscillations an artificial smoothing term can be 
added to the right-hand side of equations (10) and (11) as follows: 


predictor step 
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corrector step 
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where C is an adjustable constant. In regions of smooth flow these terms 
will be ^negligible and will not influence the solution. In regions of point- 
wise oscillations these terms will provide the effect of solution smoothing. 


Brailovskaya 

The two-step finite-difference scheme introduced by Brailovskaya in 1965 
(ref. 10) is second-order accurate in space and first-order accurate in time. 
Using this technique to difference equation (1) yields 
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The viscous terms in the predictor step are identical with the viscous terms 
in the corrector step and, therefore, need to be computed only once per time 
step. This feature reduces the required amount of computer time. The sta- 
bility requirement for the Brailovskaya scheme is the usual CFL condition 
(eq. (5)). An additional viscous stability condition is required and is given 
by equation (12). The artificial smoothing applied to the MacCormack scheme 
(eqs. (13) and (14)) was also applied to the Brailovskaya scheme. 


Boundary Conditions 

The boundary conditions described in this section were used for each 
numerical method. Three boundary conditions at both the inflow and outflow 
boundaries must be specified. At the subsonic inflow total pressure and total 
temperature were specified and held fixed. The third inflow boundary condition 
was obtained by requiring a zero gradient on static pressure. 
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At the outflow boundary (i = N) for the supersonic case the boundary 
conditions were 


„n 

Pn-1 ’ 


n 

^N-1 ’ ^N-1 


At the outflow boundary for the normal shock wave case the flow is sub- 
sonic; consequently, the boundary conditions are modified as follows 


u" = u" 


" ^exit, “N " “N-1 ’ Pn “ PN-1 
where p is specified and held fixed. Obtaining accurate results with such 

G XT L 

simple boundary conditions is made possible by adding constant area duct seg- 
ments at the inflow and outflow stations of the nozzle. 


These boundary conditions when applied to the modified Du Fort-Frankel 
code resulted in unstable oscillations at the boundaries. These oscillations 
were eliminated by two different methods. The first method was to apply 
second-order damping given by 


f ^ ^ - 2U" + U" 1 

'i 16 \ ^i+1 ^^i ^ ^i-i; 


at i = 2 for the inflow and i = N -1 for the outflow. 

The second method of removing the oscillations consisted of replacing the 
original boundary conditions with a new set given by 
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where equations (20) and (21) were used for the no shock wave case and 
equations (20) and (22) for the normal shock wave case. For shock-free flow 
solutions the second method of removing the oscillations produced the best 
results and are presented in the next section. For the normal shock wave 
case, both methods produced similar results. 
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DISCUSSION OF RESULTS 


Shock-Free Solution 


The initial condition solution for the isentropic calculation was 
established by first computing the inflow and outflow endpoints from one- 
dimensional isentropic theory. Then linear distributions for all the flow 
variables were computed between the endpoints. 

Table I sutmiarizes the computing statistics of all the results presented. 
The results of the isentropic (shock free) calculation are presented in figure 
1. Included are Mach number, pressure, and temperature distributions along 
the nozzle axis for all four numerical techniques. Overall the agreement is 
very good. In particular, all four numerically predicted values of pressure 
at the throat lie within 0.5 percent of the theoretical value. The largest 
disagreement occurs at the outflow where theory predicts an exit Mach number 
of 1.925, The numerically predicted exit Mach numbers are below this and lie 
between 1.907 and 1.920. 


The maximum error (ERR) versus the central processor unit time (CPU time) 
is presented in figure 2 where 
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The CPU time required for computing initial conditions and solution input/ 
output has been subtracted from the CPU time displayed in figure 2. The 
curves have been continued until the maximum error dropped below 0.001 
although the actual calculations were carried to 0.0001 accuracy. For the test 
problem the modified hopscotch technique is clearly the fastest of the four 
techniques tested, being 2.2 times faster than Brailovskaya, 2.5 times faster 
than MacCormack, and 4.0 times faster than modified Du Fort-Frankel . 


Normal Shock Wave Solution 

The initial conditions for the cases with a standing normal shock wave 
were obtained from one-dimensional isentropic theory. The initial solution 
was entirely subsonic with the standard expansion in the converging portion in 
the nozzle, the sonic condition at the throat, and subsonic compression in the 
diverging portion of the nozzle. This condition was chosen because the use of 
initial conditions with supersonic outflows caused difficulties when the out- 
flow pressure was specified. 

Mach number and pressure distributions are presented in figures 3 and 4. 

A standing normal shock wave with a pressure ratio of approximately 3.7 has 
been captured by all four methods at i == 38. The shock wave is spread over 
two to three grid points with minimal overshoots and no undershoots. The 
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Reynolds number was not low enough for any significant viscous effects to 
appear. The grid Reynolds numbers (pu Ax/u) were between 10 and 20. In 
general the agreement is quite good between the four techniques. 

Both the MacCormack and Brailovskaya results were computed with artifi- 
cial smoothing added; however, the smoothing was not required for a stable 
solution. Instead, it was used to improve the characteristics of the captured 
shock by reducing the overshoot and undershoot oscillations. 

Again the modified hopscotch technique is the fastest of the four methods 
tested (see table I, case 2) being 1.7 times faster than Brailovskaya, 1.8 
times faster than MacCormack, and 3.9 times faster than modified Du Fort-Frankel . 


Reynolds Number Effects 

Three of the techniques (modified hopscotch, MacCormack, and Brailovskaya) 
have viscous stability conditions and therefore, should exhibit smaller time 
steps and longer CPU times for lower Reynolds numbers. Two test cases were 
computed with lower Reynolds numbers by decreasing the total pressure (see 
table I, cases 3 and 4). Mach number distributions for these two cases are 
shown in figures 5 and 6. The effect of the reduced Reynolds number is 
clearly evident. In figure 5 (R = 11345 m-1 ) the shock wave is spread across 
five to six grid points, and in figure 6 (R = 2269 m"l ) the solution is so 
smeared by the physical viscosity that a shock wave cannot be recognized. The 
,grid Reynolds numbers are between 2 and 5 for case 3 and between 0.5 and 1.0 
for case 4. 

The shock position predicted by the modified Du Fort-Frankel technique 
for case 3 (see fig. 5) is in slight disagreement with the shock position pre- 
dicted by the other three methods. This is due to the different outflow 
boundary conditions used by the modified Du Fort-Frankel technique (see eqs. 
(19) - (22)). The effect is to alter the value of exit pressure and thus, 
change the shock position. The modified Du Fort-Frankel scheme failed to 
converge for case 4. 

As expected, the viscous stability condition was more restrictive and 
therefore dominated the low Reynolds number calculations, especially case 4. 
Modified hopscotch seemed to have a slightly more severe viscous stability 
condition than MacCormack or Brailovskaya, but, possibly due to the added 
physical viscosity, actually reached a converged solution sooner in physical 
time. For instance, in case 3, modified hopscotch was 2.8 times faster than 
Brailovskaya, 3.2 times faster than MacCormack, and 6.2 times faster than 
modified Du Fort-Frankel. 

The lack of a viscous stability limit for modified Du Fort-Frankel could 
not be fully tested due to its failure to converge for case 4. The reduced 
time step ratio exhibited by modified Du Fort-Frankel for all cases is due to 
the artificial dissipation which must be added for stable operation. Hence, 
even if modified Du Fort-Frankel is not restricted by a viscous stability 
condition, it must pay the price of a reduced time step for another reason, 
regardless of Reynolds number. 
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Artificial Smoothing 


Artificial smoothing has been used in this study on three of the four 
methods tested (MacCormack, Brailovskaya, and Du Fort- Fran kel ) . To investi- 
gate the effect of smoothing, a series of Mach number distributions for three 
different values of (smoothing constant) are presented in figure 7. The 
three curves correspond to no smoothing (C^ = 0.0), moderate smoothing (Cx = 
0.2), and massive smoothing (Cx =1.0). All three curves were computed by 
the same numerical technique (Brailovskaya) and at the same flow conditions 
(R = 45374 m-1 and Pexit^Pt " Enlargements of the Mach number profiles 

around the standing normal shock wave are presented in figure 7. The no 
smoothing case spreads the shock wave across three grid points and exhibits 
pre-shock oscillations. The moderate smoothing case, likewise, spreads the 
shock over three grid points, almost identically matching the no smoothing 
shock, but without pre-shock oscillations. The massive smoothing case spreads 
the shock over four or five grid points and effectively causes a position 
shift in the shock wave. All three profiles away from the shock are in good 
agreement regardless of how much smoothing is applied. Therefore, it is 
clear that artificial smoothing in a limited amount has helped the quality 
of the solution. 


CONCLUDING REMARKS 


The modified hopscotch technique was superior in speed for all cases 
tested, being 1.7 to 2.8 times faster than the Brailovskaya technique, 1.8 to 
3.2 times faster than the MacCormack technique, and 3.9 to 6.2 times faster 
than the modified Du Fort-Frankel technique. 

All methods tested were comparable in accuracy for the cases tested, 
with or without shock waves. 

The modified hopscotch scheme seemed to have a slightly more severe 
viscous stability condition than the MacCormack or Brailovskaya schemes. 
However, for the viscous stability restricted cases, solutions computed by 
the modified hopscotch technique actually reached steady state sooner in 
physical time than any of the other techniques tested. 
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Table I.- Sunnnary of results. 


Case 1 

Isentropic Supersonic 
R =45374 


Case 2 

Normal Shock 
R =45374 
Pexit/Pt'- ’’ 


Case 3 

Normal Shock 
R =11345 
Pexit/Pf-7 


Case 4 

Normal Shock 
R =2269 


Pexit/Pt" 


^ Physical CPU ^ Physical cPU a 1 Physical CPU m Physical CPU 

n Time Time n Time Time At„., Time Time aU., ^ Time Time 

CFL (^sec) (sec) (>*sec) (sec) T-FL (jg^) CFL 


MacCormack 1.0 345 143 5.8 . 9 338 140 5.3 . 9 546 223 8.5 . 5 896 217 14.0 


Modified 

Hopscotch 


.9 395 147 2.2 1.0 565 246 3.0 .8 514 188 2.7 .3 1193 175 6.2 


Brailovskaya 1.1 334 1 52 5.2 LI 323 161 5.2 1.0 485 222 7.5 .4 1110 221 17.6 


Modified 
Du Fort- 
Frankel 


. 5 858 177 10.7 . 5 990 216 11.8 .4 1350 243 16.9 


A Modified Hopscotch 
O Brailovskaya 

□ Modified Du Fort-Franke! (t^lO, ^=.6) 
O MacCormack 


Grid index, i 

Figure 1.- Shock-free calculation 
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Figure 4.- Pressure ratio distribution (case 2). 
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Figure 6.- Mach number distribution (case 4). 



Grid Index, i 

Figure 7.- Effect of artificial smoothing. 
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